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Linear scaling methods provide total energy, but no energy levels and canonical wavefuctions. 
From the density matrix computed through the density matrix purification methods, we propose 
an order-N (O(N)) method for calculating both the energies and wavefuctions of band edge states, 
which are important for optical properties and chemical reactions. In addition, we also develop an 
O(N) algorithm to deal with doped semiconductors based on the O(N) method for band edge states 
calculation. We illustrate the O(N) behavior of the new method by applying it to boron nitride 
(BN) nanotubes and BN nanotubes with an adsorbed hydrogen atom. The band gap of various BN 
nanotubes are investigated systematicly and the acceptor levels of BN nanotubes with an isolated 
adsorbed H atom are computed. Our methods are simple, robust, and especially suited for the 
application in self-consistent field electronic structure theory. 



I. INTRODUCTION 

Traditional electronic structure algorithms calculate 
all eigenstates associated with discrete energy levels. The 
disadvantage of this approach is that it leads to a diago- 
nalization problem that has an unfavorable cubic scaling 
in the computational effort. Linear scaling density func- 
tional or Hartree-Fock methods are an essential tool for 
the calculation of the electronic structure of large systems 
containing many atoms. 1 - The key point of the success 
of most linear scaling methods is that only the density 
matrix or localized Wannier functions which span the oc- 
cupied manifold is calculated. In these O(N) methods, 
no canonical wavefunctions or eigenvalues are available. 
However, in many cases, one may be interested in some 
eigenstates, especially the states near the Fermi level, i.e., 
band edge states. For instance, from the theory of fron- 
tier orbitals, many molecular properties are determined 
by the highest occupied molecular orbital (HOMO) and 
lowest unoccupied molecular orbital (LUMO), and fron- 
tier orbitals play an important role in chemical reactions. 
On the other hand, there are some linear scaling algo- 
rithms such as the Kim-Mauri-Galli (KMG)2 method 
need the Fermi level which can be estimated from the 
HOMO and LUMO energy. 

There are some methods which can be used to obtain 
band edge states. The most popular method for cal- 
culating states near a reference energy e re f is the folded 
spectrum method. 3 However, in this method, by squaring 
the Hamiltonian, the condition number is also squared 
and thus the difficulty of solving the equation is also 
increased. To solve this problem, Tackett et alA pre- 
sented the Jacobi-Davidson method in which the con- 
dition number and difficulty in solving for the selected 
eigensolutions is the same as the original eigenvalue equa- 
tion. Unfortunately, the implementation of the Jacobi- 
Davidson method is rather involved and its application 
is not widespread. In the field of computational mathe- 
matics, the shift-and-invert Lanczos algorithm is a well- 
known method for calculating a pair of eigenvalue and 



eigenvector near a reference energy. This method was 
used by Liang et a/£ to obtain the Fermi level in the con- 
text of linear scaling Fermi operator expansion method. 
In this method, the Lanczos method is applied to the so- 
called shift-and-invert matrix, (H — e re //) _1 , where H 
and I are the Hamiltonian and identity matrices, respec- 
tively, and e re f is the reference energy. These matrices 
are not, of course, formed explicitly. Instead, each time 
the Lanczos method requires a multiplication of a vector 
v by matrix (H — e re fl)~ l , a linear solver subroutine is 
called to solve the corresponding linear systems. If these 
linear systems are solved sufficiently accurately, the con- 
vergence of the Lanczos method is typically much faster 
compared to that when the matrix H is used in the Lanc- 
zos method. The difficulty now is that accurate numer- 
ical solution of linear systems, needed on each iteration 
of the Lanczos method, can be costly. Besides the diffi- 
culties of these methods mentioned above, when they are 
applied to get the frontier orbitals, another inconvenience 
is that a reference energy e re f must be selected. 

Here in this work, we present an alternative simple 
method to get states near gap based on linear scaling 
density matrix methods. In our method, we do not need 
the reference energy. The new O(N) method is particu- 
larly useful for calculating frontier orbitals in the frame- 
work of self-consistent field (SCF) electronic structure 
theory. Using this method, we also propose a promising 
linear scaling method which can be utilized to explore 
the energetics, defective levels, and gemoetry of doped 
semiconductors . 

This paper is organized as follows: In Sec. [Til we 
present our new O(N) methods for calculating band 
edge states and dealing with doped semiconductors. In 
Sec. IIII1 we describe the details of the implementation 
and perform some test calculations to illustrate the Tight- 
ness, robustness, and linear-scaling behavior of our meth- 
ods. In Sec. HVi we use our new methods to calculate the 
band gap of boron nitride (BN) nanotubes and the accep- 
tor level of a single H adsorbed BN nanotubes. Finally, 
our concluding remarks are given in Sec. fVl 
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II. THEORY 
A. Calculation of band edge states 

Within our method, we must first obtain density ma- 
trix p corresponding a given Hamiltonian H before we 
proceed to calculate band edge states. However, it is 
not an inconvenience in the framework of linear scaling 
SCF electronic structure theory. In principles, any lin- 
ear scaling density matrix methods can be used to obtain 
the density matrixi&l Moreover, O(N) localized orbital 
based methods can also be used to construct the density 
matrixi^S 

In the representation of molecular canonical orbitals, 
density matrix p and Hamiltonian H are diagonal matri- 
ces of the following forms: 

p = diag(l,l, ...,1,0,0,..., 0), ,^ 
H = diag(e 1 ,e 2 ,...,e Ne / 2 ,^N e /2+i,---,^N b ), 

where N e is the number of electrons of a closed-shell sys- 
tem, and N is the number of basis functions. Without 
loss of generality, we assume that 

ei < £2 < • ■ • < £NJ2 < £NJ2+1 < ■■■<e Nb , (2) 

then ejv e /2 and e^/2+1 will be the HOMO and LUMO 
energies respectively. It can be easily seen that: 

pH = Hp = diag(e 1 ,e 2 , ■ ■ .,£jv e /2,0, ... ,0). (3) 

If £Ar c /2 > 0, then e_/v e /2 will be the largest eigenvalue 
of pH. Otherwise, we can shift the Hamiltonian H to 
H + XI (A > 0) so that A + e N j 2 > 0. Clearly, A + 
e N c /2 is the largest eigenvalue of p(H + XI). Using the 
similar argument, we can prove that if —A + 6^/2+1 < 
0, —A + e Ne / 2 +i will be the smallest eigenvalue of (I — 
p)(H — XI). We should note that the parameter A can 
be set to be a large positive value without degrading the 
efficiency of the method. In practice, we find that it is 
usually reliable by setting A to be 1 Ry. The largest 
(smallest) eigenvalue and its corresponding eigenvector 
of p(H + XI) ((I — p)(H — XI)) can be computed easily 
using the well-known O(N) Lanczos method. Up to now, 
we discuss the problem in the representation of molecular 
canonical orbitals tp. In the representation of orthogonal 
basis orbitals 4>, molecular canonical orbitals tp can be 
expressed as 

i^i = J2^ C ^ ( 4 ) 

where the coefficient matrix C is a unitary matrix. Thus 
in the representation of orthogonal basis orbitals, density 
matrix p or and Hamiltonian H or can be calculated as: 

p or = CpC + , 

H or = CHC+. { 1 

Moreover, p or H or can also be obtained through a unitary 
transformation of pH. Since the unitary transformation 



of a matrix does not change its eigenvalues, we can see 
that the above results deduced using the representation 
of molecular canonical orbitals also hold in the represen- 
tation of orthogonal basis orbitals. The procedure for 
obtaining HOMO and LUMO states are illustrated in 
Fig. Ha). 

Since many first principles codes use non-orthogonal 
atomic orbitals, here we discuss the case of non- 
orthogonal basis. A general method is transforming the 
non-orthogonal basis to orthogonal basis. We achieve 
this by transforming the atomic orbital (AO) Hamilto- 
nian matrix H ao to an orthonormal basis using H or = 
ZH ao Z T and obtaining the AO density matrix p ao using 
p ao = Z T p or Z, where the inverse factor Z = L^ 1 , and L 
is the Cholesky factor for which S = LL T . The Cholesky 
transformation has been used in severval linear scaling 
densit matrix programs. We next show how to get wave- 
function in the non-orthogonal basis. In non-orthogonal 
basis, a generalized eigenvalue problem should be solved: 

H ao 1pao = £Slp ao , (6) 

where ip ao is the wavefunction in the non-orthogonal ba- 
sis. Given the wavefunction in the orthogonal basis ip or , 
which satisfies 

Horror = ZH ao Z T ^ or = eip ori (7) 

we have 

H ao Z T ip or = eZ~ 1 Z~ t Z T ip or 

= eSZ T ^ ori (8) 

1p ao = Z T l/j or . 

We also present another method to calculate band 
edge states in non-orthogonal basis without transform- 
ing to orthogonal basis. This method is particularly 
useful when localized orbitals based O(N) algorithms 
are employed. From p a oH ao — Z T p or H or Z~ T ', one 
can easily see that p a oH ao ipao — tipao is equivalent to 
p or H or Z~ T ij} ao = eZ^ T ip ao . Thus p ao H ao has the same 
eigenvalues as p or H or . We can also prove that p ao {.H ao + 
XS) has the same eigenvalues as p or (H or + A). Thus the 
largest eigenvalue of p a o(H ao + XS) will be e(HOMO) + A. 
We should point out that p r(H or + A) is hermitian, but 
Pao{H ao + AS") is not. However, this doesn't pose any 
problem since the Lanczos algorithm can also be used to 
get the extreme eigenvalues of a non-hermitian matrix. 
We can see that the calculation of HOMO state is simple 
since only p aol H aol and S are needed. However, the cal- 
culation of the LUMO state is a different story. We can 
easily prove that (/ — p ao S)(S~ 1 H ao — XI) has the same 
eigenvalues as (/ — p r){H r — XI) and —A + e(LUMO) 
is the smallest eigenvalue of (/ — p ao S)(S^ 1 H ao — XI). 
As can be seen, to calculate the LUMO state, besides 
Pao, H ao , and S, we must also have S^ 1 or S~ 1 H ao . The 
inverse of S is usually a formidable task. Fortunately, 
Gibson et at introduced an O(N) method to calculate 

S 1 Hao 
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B. Treatment of doped semiconductors 

To our best knowledge, most linear scaling methods 
are mainly applied to semiconductors or insulators with 
an energy gap. When the system is metallic or gapless, 
these O(N) methods fail or lose of effectiveness since these 
methods rely on the sparsity of the density matrix and 
the convergence of many of these methods is determined 
by the magnitude of band gap. Partial occupation is an- 
other obstacle for many popular linear scaling methods 
due to the non-idempotence of the density matrix. Here 
we propose an O(N) method to deal with doped semi- 
conductors where dopants or defects exist. Our method 
has the similar spirit as that proposed by Raczkowski and 
Fong in that a subspace larger than the occupied space is 
used. 11 In their seminal work, the subspace optimization 
method formulated in terms of localized nonorthogonal 
orbitals was employed. However, besides two 0(N 3 ) steps 
in the Grassmann conjugate gradient (GCG) algorithm, 
an additional 0(N 3 ) step of diagonalization is needed. 
Another problem is that when the orbital localization is 
used to acheive linear scaling, local minima might oc- 
cur in the subspace optimization method, resulting in a 
stalling of GCG algorithm during the last several SCF 
steps 

In our method, we treat the valence bands using the 
density matrix method, and other defective bands are 
calculated using our O(N) method for band edge states. 
For simplicity, we consider the cases where only an elec- 
tron or hole is present in a semiconductor, as shown in 
Fig. In case of n-type doping (Fig. |2Ja)), the total 
density matrix p can be calculated as 

P = pval + 0.5\lpN+l)(lpN+l\, (9) 

where p va i is density matrix corresponding to the valence 
band. In case of p-type doping (Fig. HJb)), the total 
density matrix p can be calculated as 

P = Pvai - 0.5|V'jv)(V , Jv|- (10) 

Both ipN+i and ipN are computed through the newly de- 
veloped O(N) method for band edge states. Using the 
block Lanczos algorithm, our method can also be used 
when several doped levels are present. In this case, the 
Fermi distribution can be used to get the occupation of 
doped levels. Since the valence band are well separated 
from the conduction band, p va i is sparse, and the calcu- 
lation of p va i can be carried out using traditional O(N) 
methods, such as the trace-correcting density matrix pu- 
rification (TC2) method^ Since canonical orbitals tpN+i 
and 4>n are usually delocalized, the total density matrix 
p is much denser than p va i- It is difficult to deal with 
the full density matrix. However, we notice that in fact 
only a small part of the full density matrix is used in 
the construction of the new Hamiltonian. Thus in prac- 
tice, we only construct a small part of the full density 
matrix. To make our O(N) method for the treatment 
of doped semiconductors more clear, we show the flow- 
chart of a typical calculation in Fig. [ljb). Our method 



is very simple and applicable to many doped systems. 
We should mention that our method is not a black-box 
method since some knowledge of the studied system must 
be known prior. For instance, we should know the doping 
type and number of doping levels. Typically, we can get 
this information from intuition or deduction from other 
smaller systems with similar characters. 



III. IMPLEMENTATION AND RESULTS 
A. Implementation 

Our newly developed method has been implemented 
in SIESTA^ a standard Kohn-Sham density-functional 
program using norm-conserving pseudopotentials and nu- 
merical atomic orbitals as basis sets. In SIESTA, periodic 
boundary conditions are employed to simulate both iso- 
lated and periodic systems. Here we use the O(N) TC2 
methods' to get the density matrix since it is very simple, 
robust, and efficient. The details about the implementa- 
tion of the TC2 method can be found in Ref. 13. 

In our O(N) method for doped semiconductors, to 
obtain atomic forces, it is necessary to get the energy 
weighted density matrix E when using atomic basis sets. 
Take the case as shown in Fig. [2ja) as an example, 

E = E va i + 0.5ejv+i|V'iV+i)(V , JV+i|, (H) 

where E va i is calculated from p va \. For energy weighted 
density matrix E, we also compute and save only a part of 
the full matrix. To speed up the calculation, we adopt the 
block Lanczos method to calculate the defect levels, since 
the vectors produced by the previous SCF step can be 
reused in the subsequent step. Usually, in the last several 
SCF steps, we don't need any matrix-vector multiplica- 
tions in the calculation of band edge states. Thus, when 
a geometry optimization is performed, the extra amount 
for computing defect levels using our O(N) method is al- 
most negligible. This contrasts to the method proposed 
by Raczkowski and Fongi&2*ii 

B. Validity and performance of the O(N) method 
for band edge states calculation 

All our calculations reported in this work are done in 
the local density approximation (LDA)^ Unless oth- 
erwise stated, the double-^ plus polarization functions 
(DZP) basis set is used in the calculations. 

We first validate our method by computing the HOMO 
and LUMO of H 2 molecule. The energies of HOMO and 
LUMO are -7.532 (-7.53257) and -1.375 (-1.37292) 
eV, respectively (values in parenthesis are results from 
the diagonalization method). We also compare the 
HOMO and LUMO wavefunctions with those from the 
diagonalization method, and find that the agreement is 
remarkable. 
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To check the performance of our method, we calcu- 
late the HOMO and LUMO of BN(5,5) nanotubes with 
different number of atoms in the supercells. The CPU 
time used is shown in Fig. [3] We can clearly see the lin- 
ear scaling behavior of our new method for both single-^ 
(SZ) and DZP basis sets. 

For the purpose of comparison, we also calculate the 
LUMO of BN(5,5) nanotube with 400 atoms using the 
folded spectrum method. The SZ basis is adopted. Since 
the performance of the folded spectrum method is very 
sensitive to the choice of the reference energy, several dif- 
ferent reference energies varying from the midgap posi- 
tion to the LUMO energy are chosen. The precision of the 
calculation is within 3 meV with respect to the value from 
the diagonalization. As shown in Fig. the CPU time 
used is very large, especially when the reference energy 
is close to the LUMO energy (the HOMO and LUMO 
energies are -7.075 and -2.577 eV respectively in the cur- 
rent computing parameters setting). Even when the the 
reference energy is chosen to be optimal, the folded spec- 
trum method is still slower by seventeen times than our 
new O(N) method (387 s v.s. 22 s). 

C. Validity and performance of the O(N) method 
for doped semiconductors 

We will take BN(8,0) zig-zag nanotubes with an ad- 
sorbed H atom as an example to illustrate the correct- 
ness and efficiency of our new method. As shown by Wu 
et al., a H atom prefers to adsorb on a B atom, and the 
system is a p-type semiconductor^ For a BN(8,0) nan- 
otube (128 atoms in the supercell) with an adsorbed H 
atom, the energy difference between our result and that 
from the diagonalization method is only 6 meV. And the 
force differences between our result and that from the 
diagonalization method do not exceed 0.6 me V/A. Both 
the energy and forces agreement validates our new O(N) 
method for doped semiconductors. We also deal suc- 
cessfully with a BN(8,0) nanotube with two adsorbed H 
atoms on two B sites, indicating that our method also 
works in case of systems with multi defect levels. 

Here we show in Fig. [5] the CPU time used in an 
ion step for supercells with different number of atoms. 
Clearly, our new method for doped semiconductors dis- 
plays a linear scaling behavior. 

IV. APPLICATIONS 

A. Band gap of BN chiral nanotubes 

Previous study showed that for small zigzag (chiral an- 
gle a = 0°) nanotubes the energy gap decreases rapidly 
with the decrease of radius, while armchair nanotubes 
(chiral angle a = 30°) almost have a constant energy 
gap^ Although previous experiments 17 indicated a pref- 
erence for zig-zag and near zig-zag BN tubes and a plau- 



sible explanation^ was proposed, a very recent high- 
resolution electron diffraction study on BN nanotubes 
grown in a carbon-free chemical vapor deposition pro- 
cess revealed a dispersion of the chiral angles J£ Thus a 
thorough knowledge of the dependence of the band gap 
upon the chirality of BN nanotubes is desirable. Chiral 
BN nanotubes usually contain large number of atoms in 
a unit cell, e.g., a BN(14,1) nanotube has 844 atoms in 
the unit cell. These nanotubes are difficult to be treated 
using traditional methods. Here we calculate systcmat- 
icly the band gap of BN nanotubes including chiral BN 
nanotubes. Whenever the system is large enough to be 
sampled using r-point, we use the new O(N) method 
for calculating band edge states. The results are shown 
in Fig. Two general trends are observed: first for 
BN nanotubes with similar radius, BN nanotubes with 
larger chiral angles have larger band gaps, secondly, for 
BN nanotubes with chiral angles close to zero, BN nan- 
otubes with larger radius have larger band gaps. In 
addition, we can see that for BN(n,m) nanotubes with 
n + m = k, the band gap of BN(n,fc — n) does not depend 
monotonously on the n value due to the competition of 
the two trends mentioned above, however, the band gap 
of BN(fc — [f ],[§]) (Here [|] denotes the maximal integer 
no larger than |) nanotube is the largest, and BN(fc,0) 
nanotube usually has the smallest band gap except that 
the band gap of BN(8,2) nanotube is small than that of 
BN(10,0) nanotube. The band gaps of some BN nan- 
otubes were reported previously and the results are in 
accord with ours , 16 ^ 19 and a more complete picture for 
the trend of the band gap of BN nanotubes is presented 
here. 



B. Acceptor level of H adsorbed BN nanotubes 

Wu et aL— investigated the adsorption of a hydrogen 
atom on zigzag BN(8,0) nanotube using a supercell con- 
taining 32 boron and 32 nitrogen atoms and found H 
prefers to adsorb on the boron atom which introduces 
an acceptor state in the gap. They also showed that the 
dispersion of the defect band is as large as 0.2 eV. Our 
test calculations in the T-only approximation also show 
that the acceptor levels of a single H adsorbed BN(8,0) 
nanotube using a 64- atoms or 128-atoms supercell are 
1.064 eV and 1.180 eV, respectively (Here, the accep- 
tor level is defined as the energy difference between the 
acceptor state and the top of the valence band) . In addi- 
tion, the adsorption energy of the H atom also depends 
on the chosen supercell: For instance, the adsorption en- 
ergy is —0.353 (—0.246) eV when using a 64-atoms (320- 
atoms) BN(8,0) supercell and the diagonalization (our 
linear scaling) method. All these facts suggest that larger 
supercells should be used to predict the properties of BN 
nanotubes with an isolated adsorbed H atom. Here with 
the O(N) method for doped semiconductors developed 
in this paper, we can treat much larger radius BN nan- 
otubes with truely isolated adsorbed H atom through us- 
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ing huge supercells. Three BN nanotubes are considered: 
BN(8,0) nanotube simulated using a supercell with 320 
atoms, BN(15,0) nanotube simulated using a supercell 
with 720 atoms, and BN(13,2) nanotube with 796 atoms 
in the unit cell. Here we show the distribution of the ac- 
ceptor state and the highest orbital of the valence band in 
Fig. \7\ Clearly, the acceptor state is a relatively localized 
state around the adsorbed H atom, which agrees with 
the result reported by Wu et al^. However, the highest 
orbital of the valence band is delocalized and mainly con- 
tributed by N 2p z orbitals. As can be seen from Fig. [SI 
BN(15,0) nanotube and BN(13,2) nanotube have similar 
radius but different chirality, and the radius of BN(8,0) 
nanotube is smaller. The calculated acceptor levels in- 
troduced by an isolated H atom are 1.184 eV, 1.557 eV 
and 1.563 eV for BN(8,0), BN(15,0) and BN(13,2) nan- 
otubes, respectively. Thus the position of the defect level 
is closer to the top of valence bands for smaller radius BN 
nanotubes, but does not depend significantly on the chi- 
rality. 



electronic structure methods. Based on the O(N) method 
for calculating band edge states, we further develop an 
O(N) algorithm to deal with doped semiconductors. In 
our methods, no reference energy is needed to obtain 
the band edge states, and they are especially suited for 
the application in SCF electronic structure theory. The 
O(N) behavior of the new methods is demonstrated by 
applying it to bare and H adsorbed BN nanotubes. The 
band gap of various BN nanotubes are investigated sys- 
tematicly and the acceptor levels of BN nanotubes with 
an isolated adsorbed H atom are calculated. Our al- 
gorithms could be generalized straightforwardly to spin- 
unrestricted systems^ such as magnetic semiconductors 
and diluted magnetic semiconductors^ 



V. CONCLUSIONS 

We present a simple O(N) method for calculating band 
edge states using the density matrix obtained from O(N) 
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FIG. 1: Schematic illustration of (a) the O(N) method for 
calculating the HOMO (the program flow for the LUMO cal- 
culation is similar except for some modifications as described 
in the text) and (b) the O(N) method for dealing with doped 
semiconductors. Here "density matrix" is abbreviated to 
"DM". 



FIG. 2: (Color online) Schematic illustration of the electronic 
structure of doped semiconductors: (a) n-typ doping and (b) 
p-type doping. 



FIG. 3: Total CPU time for calculating HOMO and LUMO 
of BN(5,5) nanotubes using the linear scaling method. Here, 
both SZ and DZP basis sets are used. All calculations were 
carried out on a 1.5 GHz Itanium 2 CPU workstation running 
RedHat Linux Advanced Server V2.1. 



FIG. 4: Total CPU time for calculating the LUMO of BN(5,5) 
nanotube with 400 atoms using the folded spectrum method 
with different reference energy. Here, SZ basis set is used. 



FIG. 5: Total CPU time for calculating of BN(8,0) nanotube 
with a H atom adsorbed on a boron atom using the O(N) 
method for doped semiconductors. Here, double-^ (DZ) basis 
set is used. 



FIG. 6: Band gap for various BN(n,m) nanotubes. BN(n,m) 
nanotubes with n + m = k are connected with a line to guide 
the eyes. 



FIG. 7: (Color online) (a) The acceptor state and (b) the 
highest orbital of the valence band of a BN(13,2) nanotube 
with an isolated adsorbed H atom. The insets show the en- 
larged plots around the adsorbed H atom. 
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